







## Load packages
library(sandwich)
library(list)
library(foreign)
library(plotly)
library(ggplot2)
library(MASS)

require(Zelig)
library(ggeffects)
library(jtools)
library(devtools) 



## set default working directory
setwd("F:/Dropbox/Research/PSRM/replication file/data")



sink("R_results.txt",split = TRUE, append=F)

trust_data <- read.dta("CHFS_2015.dta")



########################################################
##Figure 1
########################################################


base.est <- ictreg(trustN ~ 1, data =trust_data,treat="treat",J=4, method = "lm")
summary(base.est)
nrow(base.est$data)
base.pred.est<- predict(base.est, se.fit =T, avg=T, interval="confidence")
base.pred.est
base_pred_est <-base.pred.est$fit





########################################################
##Figure 4
########################################################


data.rural <- subset(trust_data, rural == 1)




list3.nrpp<- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural, data =data.rural,treat="treat",J=4,method = "ml")
summary(list3.nrpp)
nrow(list3.nrpp$data) 



########################################################
##Figure A6
########################################################


####Health Insurance
#ml
ml.health <- ictreg(trustN ~healthins +  logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml")
summary(ml.health)
nrow(ml.health$data)

#top ccoding error
top.coded.error.health <- ictreg(trustN ~ healthins + logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",error = "topcode")
summary(top.coded.error.health)

#uniform error
uniform.error.health <- ictreg(trustN ~ healthins + logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",error = "uniform")
summary(uniform.error.health)


#robust ml
robust.ml.health <- ictreg(trustN ~healthins +  logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",robust = TRUE)
summary(robust.ml.health)


#robust nls
robust.nls.health <- ictreg(trustN ~ healthins + logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "nls",robust = TRUE)
summary(robust.nls.health)



####pension insurance
#ml
ml.pension <- ictreg(trustN ~ pension_benefit+ logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml")
summary(ml.pension)
nrow(ml.pension$data)

#top ccoding error
top.coded.error.pension <- ictreg(trustN ~ pension_benefit+logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",error = "topcode")
summary(top.coded.error.pension)

#uniform error
uniform.error.pension <- ictreg(trustN ~ pension_benefit+logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",error = "uniform")
summary(uniform.error.pension)

#robust ml
robust.ml.pension <- ictreg(trustN ~pension_benefit+ logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "ml",robust = TRUE)
summary(robust.ml.pension)

#robust nls
robust.nls.pension <- ictreg(trustN ~ pension_benefit+logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat", J=4,method = "nls",robust = TRUE)
summary(robust.nls.pension)



####NPRS

data.rural <- subset(trust_data, rural == 1)


#ml
ml.nrps <- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =data.rural,treat="treat", J=4,method = "ml")
summary(ml.nrps)
nrow(ml.nrps$data)

#top ccoding error
top.coded.error.nrps <- ictreg(trustN ~ NRPP +logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =data.rural,treat="treat", J=4,method = "ml",error = "topcode")
summary(top.coded.error.nrps)

#uniform error
uniform.error.nrps <- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =data.rural,treat="treat", J=4,method = "ml",error = "uniform")
summary(uniform.error.nrps)

#robust ml
robust.ml.nrps <- ictreg(trustN ~ NRPP +logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =data.rural,treat="treat", J=4,method = "ml",robust = TRUE)
summary(robust.ml.nrps)

#robust nls
robust.nls.nrps <- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =data.rural,treat="treat", J=4,method = "nls",robust = TRUE)
summary(robust.nls.nrps)





########################################################
##Figure A12
########################################################


install_github('ChristopherLucas/MatchingFrontier@148a33bce48e442160c873ea4e0528df54964dac')
library(MatchingFrontier)


## set default working directory
setwd("F:/Dropbox/Research/PSRM/replication file/data")



## load data
## CFPS 2014, 74 cities
trust_data_cfps <- read.dta("CFPS_2014_matching.dta")


##NRPP
data.new <-na.omit(trust_data_cfps)
data.rural <- subset(data.new, rural == 1)



match.on.nrpp <- colnames(data.rural)[!(colnames(data.rural) %in% c('trust_cadre_dummy',
                                                                    'NRPP'))]


match.on.nrpp <- c('logfincome','age','male','school_yr','hukou','ccp','log_pop_10_13','log_gdp_pc_10_13','pro_rural')

form.nrpp <- as.formula(trust_cadre_dummy ~ NRPP+  logfincome+
                          age + male +school_yr+ hukou+ ccp +
                          log_pop_10_13+  log_gdp_pc_10_13+ pro_rural)

frontier.nrpp<-makeFrontier(dataset=data.rural,treatment ='NRPP',outcome='trust_cadre_dummy',match.on=match.on.nrpp)


estimates.nrpp <- estimateEffects(frontier.nrpp, 'trust_cadre_dummy ~ NRPP',
                                  mod.dependence.formula = form.nrpp,
                                  continuous.vars = c('school_yr', 'age', 'logfincome', 'log_pop_10_13','log_gdp_pc_10_13','pro_rural'),
                                  prop.estimated = .1,
                                  means.as.cutpoints = TRUE)




pdf("fig_a12.pdf", width = 7, height = 5)

plotEstimates(estimates.nrpp,
              ylab = "Estimate",
              xlab = "Number of Observations Pruned",
              cex.lab = 1.5,
              cex.axis = 1.5,
              las = 1,
              mod.dependence.col = "gray75",
              mod.dependence.border.col = "gray75",
              line.col = "black",
              xlim = c(0, 8200),
              ylim = c(-0.2, 0.2)
)
abline(0,0, col = 'black', lty = 2)


dev.off()








########################################################
##Table A5
########################################################



## load data
## data 1: CHFS 2015, 74 cities
trust_data <- read.dta("CHFS_2015.dta")



##health insurance

list2.base.health_ins <- ictreg(trustN ~ healthins + logtotal_income1+age + male+ school_yr +hukou+rural+ccp, data =trust_data,treat="treat",J=4,method = "ml")
summary(list2.base.health_ins)
nrow(list2.base.health_ins$data)

list3.base.health_ins <- ictreg(trustN ~ healthins + logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural, data =trust_data,treat="treat",J=4,method = "ml")
summary(list3.base.health_ins)
nrow(list3.base.health_ins$data)

##pension insurance

list2.base.pension_ins<- ictreg(trustN ~pension_benefit +logtotal_income1+age + male+ school_yr +hukou+rural+ccp, data =trust_data,treat="treat",J=4,method = "ml")
summary(list2.base.pension_ins)
nrow(list2.base.pension_ins$data)

list3.base.pension_ins<- ictreg(trustN ~pension_benefit +logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural, data =trust_data,treat="treat",J=4,method = "ml")
summary(list3.base.pension_ins)
nrow(list3.base.pension_ins$data)


##NRPS
data.rural <- subset(trust_data, rural == 1)


list2.nrpp<- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp, data =data.rural,treat="treat",J=4,method = "ml")
summary(list2.nrpp)
nrow(list2.nrpp$data) 

list3.nrpp<- ictreg(trustN ~NRPP + logtotal_income1+age + male+ school_yr +hukou+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural, data =data.rural,treat="treat",J=4,method = "ml")
summary(list3.nrpp)
nrow(list3.nrpp$data) 




########################################################
##Table A9
########################################################



list2.base.corrupt<- ictreg(trustN ~ logcorrup_num +  logtotal_income1+age + male+ school_yr +hukou+rural+ccp, data =trust_data,treat="treat",J=4,method = "ml")
summary(list2.base.corrupt)
nrow(list2.base.corrupt$data)

list3.base.corrupt<- ictreg(trustN ~ logcorrup_num +  logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural, data =trust_data,treat="treat",J=4,method = "ml")
summary(list3.base.corrupt)
nrow(list3.base.corrupt$data)




########################################################
##Table A10
########################################################


list2.base.income <- ictreg(trustN ~ logtotal_income1+age + male+ school_yr +hukou+rural+ccp,data =trust_data,treat="treat",J=4,method = "ml")
summary(list2.base.income)
nrow(list2.base.income$data)

list3.base.income <- ictreg(trustN ~ logtotal_income1+age + male+ school_yr +hukou+rural+ccp+ log_gdp_pc_10_13+log_pop_10_13+pro_rural,data =trust_data,treat="treat",J=4,method = "ml")
summary(list3.base.income)
nrow(list3.base.income$data)





########################################################
##Table A11
########################################################

rm(list=ls())
library("readxl")
library(list)
library(tidyverse)
library(stargazer)
library(Matching)
library(ebal)
library(xtable)
library(MASS)
library(dplyr)
library(tibble)
library(ggplot2)
library(forcats)
library(weights)
library(Hmisc)
critical <- abs(qnorm(0.025))




## set default working directory
setwd("F:/Dropbox/Research/PSRM/replication file/data")



trust_dat_cy_p = read.csv("trust_dat_cy_p.csv")


#-------------------------------------------------------------------------------------------------
# make the overall trust variables in "direct" dataset binary (b for binary)
trust_dat_cy_p$centraltrust_all_b <- ifelse(trust_dat_cy_p$centraltrust_all==1 | 
                                              trust_dat_cy_p$centraltrust_all==2, 0, 1)
trust_dat_cy_p$localtrust_all_b <- ifelse(trust_dat_cy_p$localtrust_all==1 | 
                                            trust_dat_cy_p$localtrust_all==2, 0, 1)

#### Subset the dataframe (for central) ####
# The subset only contains a (control) and b (central trust)
d.list.central <- trust_dat_cy_p[c("female","education","agegroup","ccpmember", "income",
                                   "life","pinterest", "self_index","treat","y",
                                   "college", "young", "high.income", "high.life", "high.pinterest",
                                   "high.self",
                                   "china",  
                                   "confucian_all",
                                   "high.china", 
                                   "high.confucian_all",
                                   "centraltrust_all", "localtrust_all")]
d.list.central <- d.list.central[!(d.list.central$treat =="c"),]

# Make a (control) = 0, b (treated) = 1
d.list.central$treat[d.list.central$treat == "a"] = "0"
d.list.central$treat[d.list.central$treat == "b"] = "1"


#-------------------------------------------------------------------------------------------------
#### Subset the dataframe (for local) ####

# The subset only contains a (control) and b (local trust)
d.list.local <- trust_dat_cy_p[c("female","education","agegroup","ccpmember", "income",
                                 "life","pinterest", "self_index","treat","y",
                                 "college", "young", "high.income", "high.life", "high.pinterest",
                                 "china", 
                                 "confucian_all",
                                 "high.china", 
                                 "high.confucian_all",
                                 "high.self","centraltrust_all", "localtrust_all")]
d.list.local <- d.list.local[!(d.list.local$treat =="b"),]

# Make a (control) = 0, c (treated) = 1
d.list.local$treat[d.list.local$treat == "a"] = "0"
d.list.local$treat[d.list.local$treat == "c"] = "1"




### list experiment
d.list.central$treat <- as.numeric(d.list.central$treat)

fit.list.central <- ictreg(y ~ female + education +  agegroup + ccpmember  + income +pinterest + life + self_index + china + confucian_all,
                           treat = "treat", J = 4 ,data = d.list.central, method = "ml")


summary(fit.list.central)


##direct question
direct.central <- lm(centraltrust_all_b ~ female + education + agegroup + ccpmember + income + pinterest + life + self_index+ china + confucian_all, 
                     data = trust_dat_cy_p)

summary(direct.central)






########################################################
##Table A12
########################################################
# Load Required Packages
library(list)
library(arm)

## set default working directory
setwd("F:/Dropbox/Research/PSRM/replication file/data")




# Load Data
data <- read.csv("replication_data_silent_victims.csv", sep=",")

###########################################################################################
# Create List Experiment on Sexual Assault During War
# a.  I won money in a lottery or competition. 
# b.  I was involved in an accident. 
# c.	I received help from a stranger.  
# d.	I was personally sexually assaulted. (Sensitive Item)

data$sexaussault <- ifelse(data$Questionnaire==2, data$J7, NA) 
data$sexaussault <- ifelse(data$Questionnaire==1, data$I7, data$J7)

###########################################################################################
# Create Variables

# Treatment Indicator
data$treatment <- ifelse(data$Questionnaire==2, 1, 0)

# Socio-Demographics
data$age <- data$B1a/10
data$female <- ifelse(data$B1b==2, 1, 0)
data$tamil <- ifelse(data$B7==2 | data$B7==3, 1, 0)
data$edu <- data$B3a

# War Experience
data$army <- ifelse(data$D2a==1, 1, 0) 
data$assist.army <- ifelse(data$D2b==1, 1, 0) 
data$displace <- ifelse(data$D7==1, 1, 0) 

# Region
data$northern <- ifelse(data$Province==9, 1, 0)
data$eastern <- ifelse(data$Province==8, 1, 0)



####direct questions

# Direct Item I
data$sexaussault_d <- ifelse(data$D4h==77, NA, data$D4h)
data$sexaussault_d <- ifelse(data$sexaussault_d==2, 0, data$sexaussault_d)
mean(data$sexaussault_d, na.rm=T)
d.1 <- lm(sexaussault_d ~ 1, data=data) 

# Direct Item II
data$sexaussault_d2 <- ifelse(data$D4i==77, NA, data$D4i)
data$sexaussault_d2 <- ifelse(data$sexaussault_d2==2, 0, data$sexaussault_d2)
mean(data$sexaussault_d2, na.rm=T)
d.2 <- lm(sexaussault_d2 ~ 1, data=data) 





# Multivariate Regressions of Indirect Measure of Sexual Violence 
# Estimates from linear regression (OLS) and binomial-logistic models (MLE)


# M2: MLE
mle.2 <- ictreg(sexaussault ~ tamil*assist.army + displace + female + age + edu + eastern, data=data, treat="treatment", J=3, method="ml")
summary(mle.2)


##OLS for M2
ols.2<- lm(sexaussault_d ~ tamil*assist.army + tamil+ assist.army+displace + female + age + edu + eastern, data = data)

summary(ols.2)


ols.3<- lm(sexaussault_d2 ~ tamil*assist.army + tamil+ assist.army+displace + female + age + edu + eastern, data = data)

summary(ols.3)




########################################################
##Table A13
########################################################



library(list)



## set default working directory
setwd("F:/Dropbox/Research/PSRM/replication file/data")

Guat = read.csv("Guatemala_data.csv")



## list,voter buying

VB1_ml<-ictreg(VBlist~c_Voted07+c_predisp4+c_recip+c_monitor3+c_female+c_age3+c_edu4+c_inc3+c_rural+c_indig, data=Guat,treat="VBtreat",J=3,method="ml")
summary(VB1_ml)

##direct, vote buying
VB1_direct<-lm(VBindv~c_Voted07+c_predisp4+c_recip+c_monitor3+c_female+c_age3+c_edu4+c_inc3+c_rural+c_indig, data=Guat)
summary(VB1_direct)



# list, Intimidation
IN1_ml<-ictreg(INlist~c_Voted07+c_predisp4+c_recip+c_monitor3+c_female+c_age3+c_edu4+c_inc3+c_rural+c_indig, data=Guat,treat="INtreat",J=3,method="ml")
summary(IN1_ml)

#direct, intimidation
in_direct<-lm(INindv~c_Voted07+c_predisp4+c_recip+c_monitor3+c_female+c_age3+c_edu4+c_inc3+c_rural+c_indig, data=Guat)
summary(in_direct)





sink()

